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Abstract 

The Meshless Local Petrov-Galerkin (MLPG) method is one of the popular meshless 
methods that has been used very successfully to solve several types of boundary value 
problems since the late nineties. In this paper, using a generalized moving least squares 
(GMLS) approximation, a new direct MLPG technique, called DMLPG, is presented. 
Following the principle of meshless methods to express everything "entirely in terms 
of nodes", the generalized MLS recovers test functional directly from values at nodes, 
without any detour via shape functions. This leads to a cheaper and even more accurate 
scheme. In particular, the complete absence of shape functions allows numerical integra- 
tions in the weak forms of the problem to be done over low-degree polynomials instead 
of complicated shape functions. Hence, the standard MLS shape function subroutines 
are not called at all. Numerical examples illustrate the superiority of the new technique 
over the classical MLPG. On the theoretical side, this paper discusses stability and con- 
vergence for the new discretizations that replace those of the standard MLPG. However, 
it does not treat stability, convergence, or error estimation for the MLPG as a whole. 
This should be taken from the literature on MLPG. 

Keywords: Meshless Methods, Moving least squares, Meshless Local Petrov Galerkin 
methods, MLPG, shape functions, diffuse derivatives 



1. Introduction 

The Moving Least Squares method (MLS) was introduced as an approximation tech- 
nique by Lancaster and Salkauskas [7], inspired by the pioneering work of Shepard [13] 
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and McLain [8, 9]. Since the numerical approximations of MLS arc based on a clus- 
ter of scattered nodes instead of interpolation on elements, many meshless methods for 
the numerical solution of differential equations were based on the MLS method in re- 
cent years. As an important example of such methods, we mention the Meshless Local 
Petrov-Galerkin (MLPG) method introduced by S.N. Atluri and his colleagues [1, 2, 3]. 
It is a truly meshless method in weak form which is based on local subdomains, rather 
than a single global domain. It requires neither domain elements nor background cells 
in either the approximation or the integration. 

In MLPG and other MLS based methods, the stiffness matrix is provided by integrat- 
ing over MLS shape functions or their derivatives. These shape functions are complicated 
and have no closed forms. To get accurate results, numerical quadrature with many inte- 
gration points is required. Thus the MLS subroutines must be called very often, leading 
to high computational costs. In contrast to this, the stiffness matrix in finite element 
methods (FEMs) is constructed by integrating over polynomial basis functions which are 
much cheaper to evaluate. This relaxes the cost of numerical integrations somewhat. For 
an account of the importance of numerical integration within meshless methods, we refer 
the reader to [4] . 

This paper avoids integration over MLS shape functions in MLPG and replaces it by 
the much cheaper integration over polynomials. It ignores shape functions completely. 

Altogether, the method is simpler, faster and more accurate than the original MLPG 
method. Wc use a generalized form of the MLS which directly approximates boundary 
conditions and local weak forms, shifting the numerical integration into the MLS itself, 
rather than into an outside loop over calls to MLS routines. We call this approach Direct 
Meshless Local Petrov-Galerkin (DMLPG) method. The convergence rate of MLPG and 
DMLPG seems to be the same, but thanks to the simplified computation, the results of 
DMLPG often are more precise than the results of MLPG. All of this is confirmed by 
numerical examples. 

2. Meshless methods 

Whatever the given problem is, meshless methods construct solutions from a trial 
space U whose functions are parametrized entirely in terms of nodes" [5]. We let these 
nodes form a set X := {xi, . . . ,xn}- Then the functions u of the linear trial space U are 
parametrizable by their values on X iff the linear functional Sxi, - ■ ■ , <5xn are linearly 
independent on U. This implies that there must be some basis ui, . . . , ?ijv of U such that 
the N X N matrix of values Uj{xk) is invertible, but we are not interested in knowing or 
constructing this basis. We only assume that the discretized problem is set up with a 
vector 

u = {u{xi),...,u{xn))'^ 
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of unknowns in "meshless" style, and all data have to be expressed in terms of these. 
Furthermore, we assume the discretized problem to consist of equations 

Afe(u) = ^fc, 1 < fc < M, (2.1) 

where we have M > N linear functionals Ai , . . . , Am and M prescribed real values 
. • . ,13m- Section 4 will describe how this is done for standard linear PDE problems, 

including the variations of the MLPG. 

The upshot of all meshless methods now is to provide good estimates Afe of all Afc 

using only values at nodes. This means that one has to find real numbers aj{Xk) with 

M 

Afe(u) = aj{Xk)u{xj} « Afe(u) for all fc, 1 < A; < M. (2.2) 

Putting the aj (Afe) into a.n MxN matrix A, one has to solve the possibly overdetermined 
linear system 

Au = b (2.3) 

with b= (/3i,...,/3m)^. 

Note that we do not mention shape functions at all. They are not needed to set up a 
linear system in terms of values at nodes. The goal just is to find good estimates for the 
target functionals Afc in terms of the values at nodes, e.g. via (2.2), to set up the matrix 
A. Note that in some cases, e.g. when the functionals Afc are derivatives at points, this 
is just a variation of a finite-difference approach. 

In a second stage, users might want to evaluate u at other places than in the nodes 
Xj. This is a problem of recovery of functions from discrete data values, and completely 
independent of PDE solving. There are various possibilities to do so, including the 
standard MLS with its shape functions, but we do not comment on these techniques 
here. 

3. Generalized moving least squares (GMLS) approximation 

Before we show how to discretize PDEs in the form (2.1), we focus on how to find good 
estimates of functional values X{u) in terms of nodal values u{xi), .... u{xj^). The clas- 
sical MLS approximates X{u) = u{x) from nodal values, minimizing a certain weighted 
discrete I2 norm. But in view of the previous section, we need more general functionals. 
Therefore we employ a generalized version of Moving Least Squares, adapted from [10]. 

Let u G C""(r2) for some to > 0, and let {iJ.j{u)}jL^ be a set of continuous linear 
functionals /i^ from the dual C"' (il)* of C"™(0). For a fixed given functional A G C™(0)*, 
our problem is the approximate recovery of the value X{u) from the values {Mj(u)}|^]^. 
This fits into the preceding section for A = Afc, 1 < A; < M and iJ.j{u) = u{xj), 1 < j < N. 
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The functionals A and /ij, 1 < j < A'^, can, for instance, describe point evaluations of 
w, its derivatives up to order m, or some local integrals that contain u or its derivatives 
in their integrands. In particular, we shall use functionals of the form (4.4) arising from 
local weak forms, or simple point evaluation functionals on the Dirichlet part of the 
boundary. 

The approximation A(u) of A(u) should be a linear function of the data l^j{u), i.e. it 
should have the form 

N 

and the coefficients aj should be linear in A. We already saw this in (2.2). As in 
the classical MLS, we assume the approximation equation (3.1) to be exact for a finite 
dimensional subspace P = span{pi,p2, ■ ■ ■ ,Pq} C C""(ri), i.e. 

AT 

aj {X)iij (p) = A(p) for all pGV. (3.2) 

As in the classical MLS, we employ the standard technique of minimizing 

1 ^ 

-^a|(A)M (3.3) 

as a function of the coefficients aj{X) subject to the linear constraints (3.2), where we use 
positive weights Wi,. . . ,Wn which later will be chosen in a specific way to localize the 
approximation, provided that A is a functional acting locally. Anyway, the weights should 

depend on the fimctionals A and In most cases, the functional A will be localized at 
some point x, and then wo shall use the standard MLS weights for evaluation at x. 

The GMLS approximation A(u) to X{u) can also be obtained as X{u) = X{p*), where 
p* GV is minimizing the weighted least-squares error functional 

(^:> ~ (P)) ^'^3 ' (^•^) 

among all p G P. This problem is independent of the functional A and can be efficiently 
applied for several functionals A for fixed functionals fij. This may simplify certain 
calculations a lot, provided that several functionals have to be estimated based on the 
same local data. Details are in [10] , including error bounds for the recovery and a proof 
that (3.1) holds for X{u) = X{p*) with the optimal solution p* of (3.4) and the optimal 
solution a* (A) of (3.3). However, in meshless methods, we need more than the single 
value X{u) = X{p*), since we finally need the solution as a vector a* (A) G with N 
values a* (A), 1 < j < A''. 
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In the special case 



X{u) = {5^D"){ti.), 



(3.5) 



the derivatives of u of order \a\ are recovered. They are called GMLS approximation 
derivatives in [10]. Some authors call them diffuse derivatives, but they not "diffuse" in 
any way. They are very good direct recoveries of the derivatives of u, but not coincident 
with the corresponding derivatives of the shape functions of the classical MLS solution. 
Our GMLS approach does not even have shape functions. Instead, derivatives are esti- 
mated directly from nodal values, avoiding the inefficient detour via classical derivatives 
of shape functions. 

Note that the use of polynomials is not mandatory, and the resulting values aj{\) 
will be independent of the chosen basis of V. However, choosing a good basis of V will 
improve stability, and the following discussion shows that V should have the property 
that \{p) should be easy to evaluate for p E V. 

Even if a different numerical method is used to minimize (3.3) or (3.4), the optimal 
solution a* (A) e can be written as 



where W is the diagonal matrix carrying the weights Wj on its diagonal, P is the N x Q 
matrix of values iij{pk), ^ < j < N, I < k < Q, and X{p) e is the vector with values 
A(pi), . . . , X{pq) of A on the basis of P. Thus it suffices to evaluate A on the space V, not 
on a certain trial space spanned by certain shape functions. This will significantly speed 
up numerical calculations, if the functional A is complicated, e.g. a numerical integration 
against a test function. Standard examples are functional of the form 



where i is a linear differential operator preserving polynomials or just the identity, and w 
is some polynomial test or weight function. Such functionals will arise for PDE problems 
in weak form in the next section. Then our generalized MLS technique will perform 
integration only over polynomials, if we use polynomials as the space V. Note that this 
generalizes to any type of functional: we finally just have to evaluate it on a polynomial. 
No other calls to MLS routines are necessary, because we do not apply the functional to 
shape functions. 

4. Problems in local weak forms 

We now write linear PDE problems in the discretized form (2.1), with special emphasis 
on the Meshless Local Petrov Galerkin Method. 



a* (A) = WP'^iPW P'^y^Xip) 



(3.6) 
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Although the technique proposed in this papcT can be used for a wide class of PDEs, 
we illustrate our approach for the Poisson problem 

Au{x) = f[x), X €Vt, 

u{x) = Unix), xGTd, (4.1) 
§^{x) = UNix), x&Tn 

where / is a given source function, the bounded domain C M"^ is enclosed by the 
boundary T = Td UTn, ud and un are the prescribed Dirichlet and Neumann data, 
respectively, on the Dirichlet boundary Tn and on the Neumann boundary r^v, while n 
is the outward normal direction. 

The simplest way of discretizing the problem in the form (2.1) is direct and global 
collocation. In addition to the trial nodes Xi, . . . , Xn for obtaining nodal solution values, 
we can choose finite point sets 

Fa CO, YnCTo, Yr, CT^, Y -YaUYnUY^, \Y\ = M 

and discretize the problem by M functionals 

Ai(w) = Au(zi) = f{zi), Zi eYn c O, 

Xj{u) = u{zj) = uoizj), Zj €YdC Tn, (4.2) 

^k{u) = §^{zk) = UN{zk), Zk&YN C Tn 

using some proper indexing scheme. In MLPG categories, this is MLPG2 [1, 2]. All 
functionals are local, and strong in the sense that they do not involve integration over 
test functions. 

For FEM style global weak discretization, one can keep the second and third part 
of (4.2), but the first can be weakened using the Divergence Theorem. With sufficiently 
smooth test functions Vi, we get 

\i{u) := I (Vu • n)vi dV — I Vu ■ Vvi dVt = I fvi dO 
Jr Jn Jn 

as a replacement of the first functionals in (4.2), leading again to (2.1). 

Following the original MLPG method, instead of transforming (4.1) into a global 
weak form, we construct weak forms over local subdomains fl^ which are small regions 
taken around nodes y in the global domain O. The local subdomains could theoretically 
be of any geometric shape and size. But for simplicity they are often taken to be balls 
B{y,a) intersected with O and centered at y with radius a, or squares in 2D or cubes in 
3D centered at y with sidelength cr, denoted by S(i/, a) D Q. The variable a parametrizes 
the local subdomain's size, and we denote the boundary within SI by F^ := Sln9Sl^. We 
call a node y internal, if the boundary of the local subdomain O^ does not intersect 
F. 
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The derivation of the local weak form starts with the local integral 

/ (Au-f)vdQ = 0, (4.3) 
Jnl 

where v is an appropriate test function on Q^. Employing the Divergence Theorem, we 
get an equation 

Xyav{u):= / (Vu-n)wrfr- / Vu-Vv dn= fv dn- unv dT (4.4) 

of the form (2.1). For nodes whose subdomain boundary does not intersect Fjv, the 
second term on the right hand side vanishes. 

Note that neither Lagrange multipliers nor penalty parameters are introduced into 
the local weak form, because the Dirichlet boundary conditions are imposed directly 
using the second line of (4.2) for suitable collocation points, usually taking a subset of 
the trial nodes. 

Some variations of MLPG differ in their choice of functionals (4.4). If the test function 
V is chosen to vanish on \ Fiv, the first integral in (4.4) is zero, and we have MLPGl. 
If the local test function v is the constant 1, the second integral vanishes, and we have 
MLPG5. 

5. Implementation 

In this section, we describe the implementation of GMLS approximations to solve the 
Poisson problem (4.1) using the weak form equations (4.4). 

At first we fix m, the maximal degree of polynomial basis functions we use. These 
form the space V := of d-variate real-valued polynomials of degree at most m. The 
dimension of this space \s Q = ('"j"'') • If the problem has enough smoothness, m will 
determine the convergence rate. 

Then we choose & set X — {xi,X2, ...,xn} C of scattered trial points which fills 
the domain reasonably well, without letting two points come too close to each other. 
To make this more precise, we need the quantities fill distance and separation distance 
which are important to measure the quality of centers and derive rates of convergence. 
For a set of points X = {xi, X2, ...,xn} in a bounded domain Cl C W^, the fill distance is 
defined to be 

hx,a = sup min \\x - Xj\\2, 
and the separation distance is defined by 

Qx = lmm\\xi- Xj\\2. 
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A set X of data sites is said to be quasi-uniform with respect to a constant Cqu > if 



In this sense, we require the set X of trial nodes to be quasi-uniform. 

We now have to define the functionals Ai , . . . , Am discretizing our PDE problem. 
This requires a selection between MLPGl, MLPG2, and MLPG5, and the decision to 
use ovcrsampling or not, i.e. M > N or M = N. Oversampling will often increase 
stability at increased cost, but wo found that in our examples it was not necessary. Since 
we have to execute the GMLS method for each functional Afe, approximating it in terms 
of function values at the trial nodes in B{yk, S) r\X, we have to make sure that the GMLS 
does not break down. This means that the matrix B of (3.6) must have rank Q, if formed 
for the nodes in B{yk,5) (1 X. In general: 

Definition 5.1. A set Z of pairwise distinct points in M.'^ is called F^-unisolvent if the 
zero polynomial is the only polynomial from which vanishes on Z. 

To give a sufficient condition for unisolvency, we need 

Definition 5.2. A set C M*^ is said to satisfy an interior cone condition if there exists 
an angle 9 G (0,7r/2) and a radius r > such that for every x G fl a, unit vector ^(a;) 
exists such that the cone 



is contained in fl. 

Theorem 5.3. ([12], see also [14]) 

For any compact domain in M'* with an interior cone condition, and any m > there 
are positive constants ho and cq such that for all trial node sets X with fill distance 
hx,n < ^0; all test points y e and all radii 6 > Cohx,n, the set B{y,5) (1 X is 
F^-unisolvent. 

This means that the placement of test nodes and the choice of weight function sup- 
ports can be linked to the fill distance of the set of trial nodes. Oversampling never 
causes problems, if the weight function support radius is kept proportional to the fill 
distance of the trial nodes. 

Some test nodes should be scattered over the Dirichlet boundary Td to impose the 
Dirichlet boundary conditions. Like in the collocation case, we denote the subset of these 
points by Yd C F n Tn- For MLPG2, we similarly define Yat, and then the setup of 
the functionals simply follows (4.2), with or without oversampling. In principle, the sets 
Y(i, Yjv, Yd need not be disjoint. 



qx < hx,Q < CquQX- 



(5.1) 



C{x, 9, r):={x + ty:y£ \\y\\2 = I, y^^ > cos 9, t G [0, r 



•]} 
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For weak problems in MLPGl or MLPG5 form, we just implement the functionals 
\k,<^k,vk of (4-4) as described in Section 4. Altogether, we follow Section 2 by imple- 
menting (2.1) via (2.2), and ending with the system (2.3). 

The order of convergence of the approximated functional to its exact value is important 
in this case. Applying the same strategy presented in [10] for Xy^ct{u) ■= D"'u{y), we can 
prove 

Theorem 5.4. Let 

A(u) = \y,„,nj{u) := [ w{x) Lu{x)dx, K = Vtl or dfll, y G O^. 

In the situation of Theorem 5.3, define Q,* to he the closure of[J^^QB{x,d). Define 

N 

where ^^(A) are functions derived from the GMLS approximation in (3.6). Then there 
exists a constant c > such that for all u € C""+^(0*) and all quasi-uniform X G ^ 
with hx,a < ho we have 

X{u) - >^)\\^^^^^ < ch^%'-lu\cr.+.^n^), (5.2) 

providing Jj^ \w{x)\dx < oo and if X{u) ^ 0, w{x) Lx^dx ^ fA(a;") 7^ Oj for some a 
with \a\ = m. Here n is the maximal order of derivatives involved in linear operator L 
and |u|cm+i(n.) := max|„|=„+i 

However, we cannot guarantee that the system (2.3) has full rank, since we only made 
sure that the rows of the system can be calculated via the GMLS if Theorem 5.3 applies. 
Oversampling will usually help if the system causes problems. 

After the solution vector u of (2.3), consisting of values u{xj) of values at the trial 
nodes is determined by solving the system, other values of the solution function u{x) 
(and also its derivatives) can be calculated in every point x G il again using the MLS 
approximation. 

Since we have direct approximations for boundary conditions and local weak forms, 
this method is called direct meshless local Petrov-Galerkin (DMLPG) method. It comes 
in the DMLPGl, DMLPG2, and DMLPG5 variations. 

In contrast to MLPG2, if the GMLS derivatives ("diffuse" derivatives) [10] are used 
instead of the standard derivatives of MLS shape functions, we have DMLPG2. As 
investigated in [10], the accuracies for calculating the matrix A of (2.3) are the same, 
but the computational cost of DMLPG2 is less. When looking into the literature, we 
found that DMLPG2 coincides with the Diffuse Approximate Method (DAM) [11]. But 



since we avoid using the word diffuse because there is nothing "diffuse" about these 
derivatives [10], we will call the method DMLPG2 or direct MLS collocation (DMLSC) 
method. 

As we saw in Section 3, in DMLPGl/5 methods the integrations are done only over 
polynomials, if the latter are used in the GMLS. For every functional Afc, 1 < k < M > N, 
the GMLS routine is called only once. There are no calls to produce vahies of shape 
functions. The standard MLPG/MLS technique at degree m implements numerical inte- 
gration by calling shape function evaluations, and thus the MLS routine is called approx- 
imatively M K times where K is the average number of integration points. Moreover, in 
standard MLPG methods the derivatives of MLS shape function must also be provided, 
while DMLPG has no shape functions at all. Consequently, DMLPG is considerably 
faster than MLPG. In addition, due to the error analysis presented in Theorem 5.4 for 
the new GMLS method, the final accuracies of both MLPG and DMLPG methods are 
the same. We will see experimentally that DMLPG is even more accurate than MLPG. 

As highlighted in [4], numerical integration in FEM is simple because the integrands 
of the elements of the stiffness matrix are polynomials. In contrast to this, the shape 
functions used in standard meshless methods are much more costly to evaluate, making 
numerical integration a much bigger challenge than for the FEM. In MLPG methods, 
numerical integrations are simpler than for various other meshless methods, since the 
local weak form breaks everything down to local well-shaped subdomains. However, 
since the integrands are based on MLS shape functions and their derivatives, a Gauss 
quadrature with many points is required to get accurate results, especially when the 
density of nodes increases. Overcoming this drawback is a major advantage of DMLPG 
methods, because the integrations are done over polynomials, like in FEM. 

It is interesting to note that if local sub-domains are chosen in DMLPG5 as S{x,a) 
(square or cube), a (d — 1) times [^] -points Gauss quadrature gives the exact solution 
for local boundary integrals around the nodes in the interior of fi. In DMLPGl, if again 
S{x, a) is chosen as a local sub-domain and if a polynomial test function is employed, 
a d-times -points Gauss quadrature is enough to get exact interior local 

domain integrals. Here, n is the degree of the polynomial test function. As a polynomial 
test function on the square or cube for DMLPGl with n = 2, we can use 



v{x;Xk) 



/ 4 \ 

i=l ^ ' 

0, otherwise 



where x = and Xk = (Cfeii •••iCfed)- In DMLPGl with balls as sub-domains, 

weight functions of the form function 

= 0(^^), (5.3) 
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with 5 = a can be used as test functions. Both of these test functions vanish on T^'XTj^, 
as required in DMLPGl. 

Note that, if the second weak forms (Green forms) are taken over local sub-domains 
and a modified fundamental solution is used as test function, the process gives the 
DMLPG4 rather than MLPG4 or the meshless LBIE method presented in [15]. In 
DMLPG4, it is better to use balls as local sub-domains, because in this case the modified 
fundamental solution, used as a test function, can be derived easily. But the test function 
is not a polynomial. 

The trial and test functions in both MLPG3 and 6 come from the same space and 
thus they are Galerkin type techniques. If we formulate the analogues DMLPG3/6, the 
integrands include shape functions again. Therefore, they annihilate the advantages of 
DMLPG methods with respect to numerical integration, and we ignore them in favour 
of keeping all benefits of DMLPG methods. 

Instead, we add some remarks about selecting m, the degree of polynomial basis 
functions in the GMLS. For m = 1, the variants DMLPG 1, 4, and 5 will necessarily 
fail. The background is that the GMLS performs an optimal recovery of a functional 
A in terms of nodal values, and the recovery is exact on a subspace V, using minimal 
coefficients at the nodal values. Thus, in all cases where the functional is zero on V 
by some reason or other, the recovery formula will be zero and will generate a zero 
row in the stiffness matrix. This happens for all variations based on functionals (4.4) 
and functionals extracted from the second weak form on interior points, since all those 
functionals are reformulations of 



and thus vanish on harmonic functions w, in particular on linear functions. Thus, for 
solving inhomogeneous problems, users should pick spaces V of non-harmonic functions, 
if they employ GMLS with exactness on V. This rules out polynomials with degree 
m < 1. 

Another closely related point arises from symmetry of subdomains. Since polynomials 
in a ball B{x,a) or a cube S{x,a) have symmetry properties, the entries of stiffness 
matrices in rows corresponding to internal points will often be the same for m = 2k and 
m = 2k + 1. Thus convergence rates often do not increase when going from m = 2 to 
m = 3, for instance. But this observation affects MLPG and DMLPG in the same way. 

6. Stability and convergence 

For the classical MLS and the generalized MLS from [10] and Theorem 5.4 it is known 
that the recovery X{u*) of values of functionals A on a true solution u* has an error of 
order 0(/i'"+^~'^), if h is the fill distance of the trial nodes, m is the degree of polynomials 
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used locally, if the exact solution u* is at least C™+^ , k is the maximal order of derivatives 
of u* involved in the functional, and if numerical integration has an even smaller error. 
In particular, the classical MLS and the new GMLS produce roughly the same stiffness 
matrices, but the GMLS has a considerably smaller computational complexity. 

However, the error committed in the approximation of the test functionals in terms 
of function values at nodes does not always carry over to the convergence rate of the 
full algorithm, since there is no stability analysis, so far. Even if perfect stability would 
hold, the best one can expect is to get the convergence rate implied by the local trial 
approximation, i.e. by local polynomials of degree m. This would again mean a rate 
of C'(/?7"+^~'^), but only if the solution is indeed locally approximated by polynomials 
of that degree. In fact, the next section will show that this rate can often be observed. 
But our symmetry arguments at the end of the previous section show that sometimes 
the degree m = 2fc + 1 cannot improve the behavior for m = 2k, because the odd-degree 
polynomials simply do not show up in most of the calculations for the stiffness matrix. 

7. Numerical results 

In this section some numerical results are presented to demonstrate the efficiency of 
DMLPG methods and its advantages over MLPG methods. We consider the Poisson 
equation (4.1) on the square [0, 1]^ C with Dirichlet boundary conditions. Since we 
want to study convergence rates without being limited by smnoothness of the solution, 
we take Franke 's function [6] 



where (x, y) denotes the two components of .x G M^, as analytical solution and calculate 
the right hand side and boundary conditions accordingly. Note that Franke's function is 
a standard test function for 2D scattered data fitting. Regular mesh distributions with 
mesh-size h are provided in all cases, though the methods would work with scattered 
data. We do not implement oversampling in the results of this paper. In fact, the trial 
and test points are chosen to be coincident. Also, the shifted scaled polynomial 



where z is a fixed evaluation point such as a test point or a Gaussian point for integration. 



[10], it is shown that this choice of basis function leads to more stable and accurate MLS 





is used instead of the natural polynomial basis {•x"}o<|Q:|<m for MLS approximation. In 
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approximation. We use the shifted scaled basis for both MLPG5 and DMLPG5 methods 
with m = 2,3 and 4. The Gaussian weight fimction 



f exp ( - (||x - x,\yc)^) exp ( - {S/cf) 



exp(-(Vc)2) 

0, l|a;-a;oll2 > (5 



is employed where c = cqH is a constant controlling the shape of the weight function and 
6 = 6oh is the size of the support domains. 

Let m = 2 and set cq = 0.6 and 6o = 2m. At first the local sub-domains are taken 
to be circles. To get the best results in MLPG we have to use an accurate quadrature 
formula. Here a 20-points regular Gauss-Legendre quadrature is employed for numerical 
integrations over local sub-domains. 

Numerical results, for different mesh-sizes h, are presented in terms of maximum 
errors, ratios and CPU times used for MLPG5 and DMLPG5 in Tables 1. 



Table 1. The 


maximum errors, 


ratios and CPU times used for MLPG5 and DMLPG5 with m 


= 2 


h 


MLPG5 




DMLPG5 




CPU time used 




II £|| OO 


ratio 


II ^11 oo 


ratio 


MLPG5 


DMLPG5 


0.2 


0.44 X IQ-i 




0.23 X 10-1 




0.4 sec. 


0.2 sec. 


0.1 


0.15 X 10-1 


1.59 


0.72 X 10-2 


1.68 


1.2 


0.3 


0.05 


0.73 X 10-2 


0.99 


0.20 X 10-2 


1.84 


6.5 


1.4 


0.025 


0.24 X 10-2 


1.61 


0.58 X 10-3 


1.80 


68.5 


6.5 


0.0125 


0.66 X 10-3 


1.85 


0.14 X 10-3 


1.98 


2016.0 


52.1 



The mesh-size h is divided by two row by row, therefore the ratios are computed by 

leWlloc 



l0g2 



||e(V2)||c 



Both methods have nearly the same order 2, which cannot be improved for this trial 
space, since the expected optimal order ism+1 — A; = 3 — 1 = 2. But significant 
differences occur in the columns with CPU times. As we stated before, this is due to 
restricting local integrations to polynomial basis functions in DMLPG rather than to 
integrate over MLS shape functions in the original MLPG. We could get the same results 
with fewer integration points for DMLPG, but to be fair in comparison, we use the same 
quadrature. 

In addition, to give more insight into the errors, the maximum errors of MLPG5 and 
DMLPG5 are illustrated in Fig. 1. Once can see that DMLPG is more accurate, maybe 
due to avoiding many computations and hence many roundoff errors. 
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Figure 1: Comparison of MLPG5 and DMLPG5 in terms of maximum errors for m = 2. 

In Table 2 and Fig. 2, we have compared MLPG5 and DMLPG5 in case m = 3. The 
convergence rate stays at 2 for both methods, since the third-degree polynomials cannot 
contribute to the weak form. The figure shows approximately the same error results. 
But the CPU times used are indeed different. 

Table 2. The maximum errors, ratios and CPU times used for MLPG5 and DMLPG5 witii m = 3 
MLPG5 DMLPG5 CPU time used 



h ||e||oo ratio ||e||oo ratio MLPG5 DMLPG5 



0.2 


0.28 X 10^1 




0.23 X 10-1 




0.8 sec. 


0.2 sec. 


0.1 


0.13 X 10-1 


1.08 


0.74 X 10-2 


1.62 


1.8 


0.4 


0.05 


0.33 X 10-2 


1.98 


0.20 X 10-2 


1.89 


9.7 


1.6 


0.025 


0.78 X 10-3 


2.09 


0.58 X 10-3 


1.80 


87.7 


7.6 


0.0125 


0.19 X 10-3 


2.06 


0.15 X 10-3 


1.98 


2293.3 


56.1 



In the results presented up to here, we used balls (circles) as local sub-domains. Now 

we turn to use squares for both MLPG5 and DMLPG5. Also, we run the programs with 

TO = 4 to see the differences. The parameters co — 0.8 and 5q = 2m are selected. In 

DMLPG5, a 2-points Gaussian quadrature is enough to get exact numerical integrations. 

But for MLPG5 and the right hand sides we use a lO-points Gaussian quadrature for 

every sides of squares. The results are depicted in Table 3 and Fig. 3. DMLPG is more 

accurate and approximately gives the full order 4 in this case. Beside, as we expected, 

the computational cost of DMLPG is remarkably less than MLPG. 
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Figure 2: Comparison of MLPG5 and DMLPG5 in terms of maximum errors for m = 3. 



Table 3. The 


maximum errors, 


ratios and CPU times used for MLPG5 and DMLPG5 witli m 


= 4 


h 


IVILPG5 




DMLPG5 




CPU time used 




II oo 


ratio 


||e||oo 


ratio 


MLPG5 


DMLPG5 


0.2 


0.10 X 10" 




0.12 X 10° 




0.5 sec. 


0.2 sec. 


0.1 


0.25 X 10~i 


2.04 


0.17 X 10-1 


2.87 


2.7 


0.2 


0.05 


0.78 X 10-2 


1.66 


0.12 X 10-2 


3.75 


19.2 


0.9 


0.025 


0.79 X 10-3 


3.30 


0.75 X 10-1 


4.04 


142.2 


4.7 


0.0125 


0.55 X 10-* 


3.86 


0.43 X 10-^ 


4.12 


2604.9 


43.9 



Results for MLPGl and DMLPGl turn out to behave similarly. As we know, MLPGl 
is more expensive than MLPG5 [1, 2], but there is no significant difference between 
computational costs of DMLPG5 and DMLPGl. Therefore the differences between CPU 
times used for MLPGl and DMLPGl are absolutely larger. 

All routines were written using Matlab® and run on a Pentium 4 PC with 2.50 GB 
of Memory and a twin-core 2.00 GHz CPU. 



8. Conclusion 

This article describes a new MLPG method, called DMLPG method, based on gener- 
alized moving least squares (GMLS) approximation for solving boundary value problems. 
The new method is considerably faster than the classical MLPG variants, because 
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Figure 3: Comparison of MLPG5 and DMLPG5 in terms of maximum errors for m = 4. 



• direct approximations of data functionals are used for Dirichlet boundary conditions 
and local weak forms, 

• local integrations are done over polynomials rather than over complicated MLS 
shape functions, 

• numerical integrations can sometimes be performed exactly. 

The convergence rate of both methods should be the same, but thanks to avoiding many 
computations and roundoff errors, and of course by treating the numerical integrations 
in a more elegant and possibly exact way, the results of DMLPG turn often out to be 
more accurate than the results of MLPG. 

On the downside, DMLPG does not work for m = 1 since it locally uses (harmonic) 
linear functions instead of complicated shape functions. But most MLPG users choose 
higher degrees anyway, in order to obtain better convergence rates. 

Altogether, we believe that the DMLPG methods have great potential to replace the 
original MLPG methods in many situations. 
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